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[i]  Relativistic  (>1  MeV)  electron  flux  increases  in  the  Earth’s  radiation  belts  are 
significantly  underestimated  by  models  that  only  include  transport  and  loss  processes, 
suggesting  that  some  additional  acceleration  process  is  required.  Here  we  use  a  new, 
three-dimensional  code  that  includes  radial  difliision  and  quasi-linear  pitch  angle  and  energy 
diffusion  due  to  chorus  waves,  including  cross  terms,  to  simulate  the  9  October  1990 
magnetic  storm.  The  diffusion  coefficients  are  activity  dependent,  and  time-dependent 
boundary  conditions  are  imposed  on  all  six  boundary  faces,  taken  from  fits  to  CRRES 
Medium  Electrons  A  electron  data.  Although  the  main  phase  dropout  is  not  fully  captured, 
the  persistent  phase  space  density  peaks  observed  during  the  recovery  phase  are  well 
explained,  but  this  requires  both  chorus  wave  acceleration  and  radial  diffusion. 

Citation:  Albert,  J.  M.,  N.  P.  Meredilh,  and  R.  B.  Home  (2009),  Three-dimensional  diffusion  simulaiion  of  outer  radialion  bell 
clecirons  during  ihe  9  Ociober  1990  magneiic  storm,  J.  Geophys.  Res.,  1 1 4,  A09214,  doi:10.1029/2009JA0l4336. 


1.  Introduction 

[2]  Outer  zone  radiation  belt  electrons  exhibit  highly 
dynamic  behavior  during  geomagnetic  storms.  It  has  been 
well  documented  that  the  energetic  flux  drops  rapidly 
during  the  storm  main  phase  but  recovers  over  several  days, 
often  to  higher  than  original  levels.  Radial  diffusion  accel¬ 
erates  particles  (as  they  move  inward  at  constant  first  and 
second  adiabatic  invariant)  but  is  hard  pressed  to  reproduce 
the  rate  and  extent  of  the  recovery,  especially  when  losses 
are  considered,  without  an  additional  source  of  energization. 

[3]  The  moderate  storm  that  occurred  on  9  October  1990 
has  been  particularly  well  studied  because  of  its  detailed 
observation  by  CRRES.  Brautigam  and  Albert  [2000] 
simulated  it  with  activity-dependent  radial  diffusion  and  a 
realistic,  variable  outer  boundary  condition.  Plasmaspheric 
hiss  was  the  only  loss  process  considered.  This  model  was 
found  to  work  reasonably  well  for  electrons  with  first 
adiabatic  invariant  M  %  100-300  MeV/G  but  was  unable 
to  account  for  the  increase,  and  inward  pointing  phase  space 
density  gradient,  for  M  ^  700-1000  MeV/G.  Many  other 
one-dimensional  simulations  of  radial  diffusion  have  been 
performed,  usually  with  time  scale  estimates  for  wave- 
induced  losses  [e.g.,  Shprits  and  Thorne,  2004;  Shprits  et 
ai,  2005,  2006b;  Fei  et  al.,  2006;  Lam  et  al.,  2007]  and 
with  an  estimated  internal  source  term  [Tm  et  al.,  2009].  The 
results  are  mixed  but  generally  support  the  finding  that 
radial  diffusion  is  insufficient. 

[4]  Energy  difliision,  caused  by  cyclotron  resonant  inter¬ 
actions  with  whistler  mode  chorus  waves,  is  frequently 
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invoked  as  a  candidate  mechanism  for  additional  energiza¬ 
tion  [Horne  and  Thorne,  1998;  Summers  et  al.,  1998]. 
Indeed,  the  gradual  aceeleration  of  electrons  to  relativistic 
energies  during  the  recovery  phase  of  the  9  October  1990 
storm  was  associated  with  prolonged  substorm  activity  as 
monitored  by  the  AE  index,  electron  injections  at  subrela- 
tivistic  energies,  and  enhanced  chorus  amplitudes  [Meredith 
et  ai,  2002].  Furthermore,  flat  topped  electron  pitch  angle 
distributions,  characteristic  of  pitch  angle  and  energy  scat¬ 
tering  by  resonant  wave-particle  interactions  with  whistler 
mode  chorus  waves,  developed  at  MeV  energies  [Horne  et 
ai,  2003].  Much  work  has  been  done  in  recent  years  to 
document  enhanced  chorus  waves  during  storms  [e.g., 
Meredith  et  al,  2003a;  Smith  et  ai,  2004],  to  evaluate  the 
corresponding  quasi-lincar  diffusion  coefficients  [Albert, 
2005;  Glauert  and  Horne,  2005],  and  to  estimate  the 
particle  evolution  using  a  onc-dimensional  (1-D)  energy 
diffusion  equation  [e.g..  Summers  and  Ma,  2000;  Summers 
et  al.,  2002;  Horne  et  ai,  2005a,  2005b]. 

[5]  Several  idealized  2-D  simulations  of  difliision  in  pitch 
angle  and  energy  near  L  -  4.5  have  been  done,  accounting 
for  chorus  waves  [Albert  and  Young,  2005;  Shprits  et  ai, 
2006a;  Tao  et  ai,  2008;  Xiao  et  ai,  2009],  hiss  combined 
with  magnetosonic  waves  [Tao  et  al.,  2009],  and  chorus 
waves  combined  with  VLF  hiss  and  electromagnetic  ion 
cyclotron  (EMIC)  waves  in  high -density  plumes  [Li  et  ai, 

2007] .  None  of  these  studies  included  radial  diffusion.  No 
studies  seem  to  have  solved  the  local  diffusion  equation 
with  radial  diffusion  treated  as  a  source  or  loss  term, 
although  Thorne  et  al.  [2007]  used  lifetimes,  obtained  from 
a  pitch  angle  diffusion  equation,  in  separate  I  -D  equations 
for  radial  diffusion  and  energy  diffusion. 

[6]  Some  preliminary  three-dimensional  simulations,  in¬ 
cluding  radial,  pitch  angle,  and  energy  diffusion,  have  been 
performed  [Varotsou  et  ai,  2005,  2008;  Subbotin  et  ai, 

2008] .  Furthermore,  some  progress  has  been  made  in  adding 
pitch  angle  and  energy  diffusion  to  advection-driven  ring 
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Figure  1.  Models  of  the  chorus  wave  magnetic  field 
intensity,  based  on  CRRES  measurements,  used  to  calculate 
pitch  angle  and  energy  diffusion  coefficients  for  different 
ranges  of  Kp. 


current  codes,  which  are  bounce  averaged  but  not  driff 
averaged  and  so  are  essentially  four  dimensional.  Fok  et 
al.  [2008]  treated  pitch  angle  and  energy  diffusion  by  chorus, 
while  Jordanova  et  al.  [2008]  included  pitch  angle  diffusion 
by  EMIC  waves. 

[7]  It  is  widely  believed  that  the  generation  of  chorus 
waves  inherently  involves  nonlinear  wave*particle  coupling 
[e.g.,  Nunn.,  1974;  Nunn  et  al.,  1997;  Katoh  and  Omura, 
2007a,  2007b].  However,  the  particles  involved  are  generally 
of  much  lower  energy  than  the  ones  considered  here,  which 
are  taken  to  interact  “parasitical ly’’  with  fully  developed 
chorus.  Individual,  coherent  whistler  mode  waves  can  lead  to 
particle  diffusion,  phase  trapping,  or  phase  bunching  (with¬ 
out  trapping),  depending  primarily  on  the  competing  effects 
of  wave  amplitude  and  background  magnetic  field  inhomo¬ 
geneity  at  resonance  [Albert,  2000,  2002;  Trakhtengerts 
et  al.,  2003;  Omura  and  Summers,  2006;  Demekhov  et  al., 
2006;  Bortnik  etai,  2008].  However,  the  applicability  of  this 
picture  to  the  global  evolution  of  energetic  particle  distribu¬ 
tions  during  storm  conditions  has  not  yet  been  established. 
This  paper  is  based  on  quasi-linear  diffusion  both  because  to 
some  extent  “it  seems  to  work”  and  as  a  basis  for  comparison 
with  future  developments  in  nonlinear  modeling. 

[s]  This  paper  combines  diffusion  by  chorus  waves  with 
radial  diffusion  in  a  carefully  chosen  three-dimensional  grid. 
Cross  diffusion,  which  was  not  treated  by  any  of  the  papers 
just  cited  except  Albert  and  Young  [2005],  Tao  et  al.  [2008, 
2009],  and  Xiao  et  al.  [2009],  is  included.  CRRES  Medium 
Electrons  A  (MEA)  data  are  used  to  evaluate  time-dependent 
boundary  conditions  at  low  and  high  E  and  no,  as  well  as  at 
high  and  low  L.  This  required  substantial  fitting,  interpola¬ 
tion,  and  extrapolation  of  the  data,  as  described  in  section  3. 
The  particle  data  were  used  to  drive  the  boundaries  only; 
affer  initialization,  data  were  not  assimilated  into  the  interior 
of  the  grid  (as  was  done  by  Shprits  et  al.  [2007]). 

2.  Diffusion  Equation 

[9]  Cyclotron  resonant  wave-particle  interactions  break 
the  first  two  adiabatic  invariants,  while  drift  resonant 


electric  and  magnetic  fluctuations  break  only  the  third 
invariant.  The  assumptions  of  continuous,  small,  uncorre- 
latcd  resonances  lead  to  a  multidimensional  diffusion  equa¬ 
tion  for  the  phase  space  density,  /,  written  as 


dt 


^  fr.  r.  ^f\  ^  fr.  \ 


(1) 


The  cyclotron  frequency  and  drift  frequency  interactions  are 
considered  uncoupled,  so  no  terms  involving  D13  or  D23  are 
included. 

[10]  It  is  common  to  change  variables  from  (Ji,  J2,  J3)  to 
(qq,  E,  L),  where  L  (often  denoted  L*)  labels  the  drift  shell 
[Roederer,  1970],  E  is  energy,  and  a©  denotes  equatorial 
pitch  angle.  Actually,  in  a  nonaxi symmetric  magnetic  field, 
the  minimum  (equatorial)  value  of  a  will  vary  for  different 
magnetic  field  lines  of  a  particle’s  drift  shell.  However,  it  is 
reasonable  to  relate  the  two  sets  of  variables  using  the 
expressions  suitable  for  a  dipole  field.  This  can  be  consid¬ 
ered  simply  a  convenient  change  of  variables,  as  long  as  the 
invariants  are  properly  computed  in  a  realistic  magnetic 
field  model,  and  gives 


Ot  GOao  V  P 


Of  ,  A 


p^  Oao 


dp)  GOp  \  p  Oao 


dL  A  dC 


(2) 


where  G  =p^T{nQ)  sinao  cos  qq  [Schulz,  1991]  and  Tfao)  ^ 
1.30-0.56  sin  ao  is  the  normalized  bounce  period  [e.g., 
Lyons  et  al.,  1972].  It  is  understood  that  the  L  derivatives 
are  evaluated  at  fixed  (J|,  not  fixed  (ao,  p).  Because  E 
and  p  are  simply  related,  terms  like  “energy  diffusion”  and 
“diffusion  in p*'  will  often  be  used  interchangeably. 

2.1.  Pitch  Angle  and  Energy  Diffusion  Coefficients 

[11]  The  pitch  angle  and  energy  diffusion  coefficients  arc 
evaluated  according  to  bounce-  and  drift-averaged  quasi- 
linear  theory  [Albert,  2005;  Glauert  and  Horne,  2005], 
requiring  models  of  wave  intensity,  B^,  and  its  distribution 
in  frequency  and  wave  normal  angle,  as  well  as  values  of 
the  plasma  frequency-to-gyrofrequency  ratio,  fpjfce-  As 
mentioned,  only  whistler  mode  chorus  waves  will  be 
considered.  Values  of  B^,  and  fpjfce  were  taken  from 
statistical  maps  of  CRRES  observations,  compiled  with 
resolution  of  1  h  in  magnetic  local  time  and  0.1  in  L.  The 
maps  were  also  parameterized  by  Kp  (into  three  ranges:  Kp 
<2,2  <Kp  <A,  and  Kp  >  4),  and  by  latitude  (“equatorial,” 
within  15®  of  the  equator,  and  “midlatitude,”  15®-30®  off 
the  equator).  The  resulting  values  are  shown  in  Figures  1 
and  2.  A  similar  model,  parameterized  by  AE,  was  presented 
by  Meredith  et  al.  [2003b]. 

[12]  The  frequency  and  wave  normal  angle  distributions 
were  represented,  as  usual,  by  truncated  Gaussians.  The 
peak,  width,  lower  cutoff,  and  upper  cutoff  for  w  and  x  =  tan 
6  were  (a;^,  8u),  ujyc)  =  (0.35,  0.15,  0.125,  0.575)0*,^ 
and  (x,„,  6x,  ^max)  =  (0,  tan  30®,  0,  1),  respectively. 
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Figure  2,  Models  of  equatorial  fpjfce>  based  on  CRRES 
measurements,  used  to  ealeulate  piteh  angle  and  energy 
diffusion  eoeffieients  for  different  ranges  of  Kp. 


These  were  used  to  eompute  tables  of  diffusion  eoeffieients 
for  89  integer  values  of  ao,  40  values  of  E  between 
0.0 1  and  10  MeV,  and  9  values  of  fpe/fer  between  I  and 
20,  using  the  eomputational  teehniques  of  Albert  [2005]. 
These  tables  were  then  sealed  in  B\.  and  interpolated  in 
jpjfccy  in  eonjunetion  with  the  statistieal  maps,  to  obtain 
drift-averaged  diffusion  eoeffieients  for  the  three  ranges 
of  Kp.  A  very  similar  proeedure  was  followed  by 
Varotsou  et  ai  [2008].  Results  at  L  -  4.55  are  shown  in 
Figure  3. 

2.2.  Radial  Diffusion  Coefficients 
[i3]  Eleetrie  and  magnetie  radial  diffusion  eoeffieients. 
Dll  =  +  D,„,  were  taken  from  Brautigam  and  Albert 

[2000]  and  are  also  Kp  dependent.  The  magnetie  eontribu- 


tion  is  given  as  =  lo®  <>.325^10 
while 


o=i£i£ _ (3) 

This  expression  is  in  Gaussian  units,  with  Bq  ^  0.31  G,  and 
is  based  on  eleetrie  field  fluetuations  with  E  and  the 
fluetuation  deeay  time  T  given  by 


£:  =  Eo  +  ^i(A>- 1)  7^  =  2700  s,  (4) 

with  numerieal  values  Eq  =  3.33  x  10'^  statvolt/em 
(O.l  mV/m)  and  Ex  =  8.67  x  10  statvolt/em  (0.26  mV/m). 
The  drift  frequeney  LOp  may  be  written  as 


3  ^  rricC  (p/mc)^  sin^  oq 

2  eBo  Rl  /,  ,  /  '/  .2  •  2 

V  •  +  sin^ao 


(5) 


Representative  values  of  Du  are  also  shown  in  Figure  3. 
2.3.  Variables  and  Grids 

[14]  Radial  diffusion  oeeurs  at  eonstant  first  and  seeond 
adiabatie  invariant,  so  it  is  most  simply  treated  in  the 
variables  (Jx,  J2,  /O.  Cyelotron  resonant  interaetions  are 
more  naturally  expressed  as  diffusion  in  piteh  angle  and 
energy,  both  beeause  the  boundaries  eorrespond  more 
elosely  to  real  partiele  deteetors  and  beeause  terms  involv¬ 
ing  eross  diffusion  are  typieally  dominated  by  piteh  angle 
diffusion.  However,  eross  difftision,  whieh  expresses  the 
physieal  relationship  between  resonant  ehanges  in  ao  and 
p,  can  have  signifieant  eonsequenees,  sinee  typieally 
D.^^njp^  >  >  Dpp-  Thus,  it  is  preferable  to  retain 

it  despite  the  numerieal  diffieulties  it  presents  to  straight¬ 
forward  finite  differeneing  in  (oo,  E,  L)  [Albert^  2004, 


10.00 
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Figure  3.  Diffusion  eoeffieients  at  L  =  4.55,  units  are  per  day.  The  sign  of  D,^^p  is  indieated  by  a  (red 
for  positive,  blue  for  negative).  Du  is  evaluated  at  Kp  =  I,  /T/?  =  3,  and  Kp  ~  5. 
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Figure  4.  The  3-D  computational  grid,  in  (a©,  E^L)  space. 


2009].  These  difficulties  may  be  overcome  in  a  number  of 
ways  [Albert  and  Young,  2005;  Tao  et  al,  2008,  2009; 
Xiao  et  ai,  2009]. 

[15]  The  method  of  Albert  and  Young  [2005]  used  the 
diffusion  coefficients  themselves  to  constructs  new  varia¬ 
bles,  (Q\,  Q2)  at  a  fixed  L,  in  which  cross  diffusion  vanished; 
it  consisted  of  choosing  Qi  =  a©  and  integrating  a  differential 
equation  for  curves  of  constant  Q2.  This  can  be  carried  out 
independently  at  each  L.  To  make  radial  diffusion  easy  to 
implement,  the  three-dimensional  grid  is  generated  from  a 
convenient  set  of  (a©,  £)  at  one  L  value,  mapped  to  other  L  at 
constant  and  Ji  as  in  a  dipole  field,  and  then  converted  to 
(Gb  Qi)  This  gives 


[16]  Since  the  points  are  not  regularly  aligned  in  the  (0i, 
02)  plane,  finite  differencing  requires  interpolation  in  02* 
though  not  in  Qx  (since  the  mapping  in  L  preserves 
alignment  in  Q\  =  a©).  Exactly  analogous  interpolation 
would  be  required  even  without  cross  diffusion,  since  the 
mapped  points  are  not  aligned  in  the  (a©,  E)  plane  either 
[Subbotin  et  al.y  2008].  This  procedure  was  carried  out 
using  diffusion  coefficients  for  each  of  the  three  Kp  ranges. 
When  the  appropriate  range  of  Kp  changes  in  the  course  of  a 
run,  the  values  of  (Gi,  Q2)  (and  Dj,  D2,  f)  are  changed,  but 
the  grid  points  retain  their  values  of  a©,  E,  L,  and / 

[17]  Figure  4  illustrates  how  mapping  in  L  shifts  and 
distorts  the  range  of  E.  AX  L  =  6.15,  the  computational  grid 
covers  E  =  0.2-2  MeV,  while  at  L  =  3.55  this  becomes 
roughly  0.5-4  MeV.  With  this  scheme,  there  are  no  wasted 
grid  points;  all  of  the  grid  points  can  couple  to  the  compu¬ 
tational  domain  through  diffusion  in  all  of  the  variables. 
Figure  5  shows  the  grid  points  in  more  detail  in  (a©,  E) 
planes  at  several  values  of  L.  Also  shown  are  the  same 
physical  points  plotted  in  (Qi,  Q2)  coordinates  (evaluated  for 
Kp  <  2),  as  well  as  in  (^1,^2)  coordinates.  Grid  points  plotted 
in  red  lie  within  the  energy  range  of  the  CRRES  MEA 
detector  so  that  actual  measurements  are  at  least  potentially 
available  to  initialize  the  flux  values  and  for  comparison 
during  the  simulation. 

3.  CRRES  MEA  Data  Processing 

[is]  The  version  of  MEA  data  used  by  Brautigam  and 
Albert  [2000]  was  limited  by  both  saturation  and  contam¬ 
ination  by  high-energy  protons,  which  prohibited  the  use  of 
the  two  lowest  energy  channels.  The  version  used  here, 
available  through  the  National  Space  Science  Data  Center, 
has  been  reprocessed,  including  a  “foldover  correction,” 
which  overcomes  these  limitations  [Vampola,  1998; 
Lemaireetai,  1998].  Thus,  this  MEA  data  set  provides 
flux  at  a  =  5,  10,  . . .,  90°,  and  17  values  of  energy 
(E  =  0.148  MeV  to  1.581  MeV),  every  60  s.  The  problem  is 
to  determine  values  at  points  (o:/,  £),  L*)  of  the  computa¬ 
tional  grid,  at  any  time  t.  (For  grid  points,  a/  means 
equatorial  pitch  angle  q©.) 

[19]  The  CRRES  ephemeris  files  used  provide  time- 
tagged  satellite  location  and  local  magnetic  field  B  every 
30  s.  For  each  entry,  the  Office  National  d’Etudes  et  de 
Recherche  Aerospatiales  code  [Boscher  et  ai,  2008]  was 
used  with  the  Olson-Pfitzer  quiet  and  International  Geo¬ 
magnetic  Reference  Field  magnetic  field  models  to  deter¬ 
mine  and  the  adiabatic  invariants  K  and  L  corresponding  to 
the  1 8  local  pitch  angle  values.  (L  depends  weakly  on  a  but 
not  on  energy.  K  is  defined  to  be  proportional  to  J2^\fT\ 
[Schulz,  1991]  so  that  it  is  also  independent  of  energy.) 
Results  with  L  within  6L  =  0.05  of  a  grid  value  were 
recorded,  along  with  the  earliest  and  latest  times,  t\  and  ^2* 
of  each  “visit”  to  each  L*.  Equatorial  pitch  angles  a©  were 
then  determined  from 


(It 
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OL  OL'  ^ 


r(sinQo)  K 
sin  tto  LR^  y/B^  ’ 


where  f  =  l^^i,  G2)|»  and  again  radial  diffusion 

operates  at  fixed  {Ji,  J2\  not  fixed  {Q\,  Qi). 


where  is  the  radius  of  the  Earth,  is  the  value  of  the 
(model)  magnetic  field  at  the  equator,  and  f(sin  a©)  is  taken 
to  be  the  function  corresponding  to  a  dipole.  (See  the 
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Figure  5.  Computational  grid  points  in  2-D  for  Kp  <  2,  expressed  in  several  sets  of  variables.  The 
red  points  are  in  the  range  of  the  CRRES  MEA  deteetor.  The  blue  eurves  indieate  several  values  of 
first  adiabatie  invariant  A/,  in  MeV/G.  The  green  eurves  indieate  the  referenee  value  J2  =  1.78  x 
10  g(em/s)/?<,. 


diseussion  in  seetion  2.)  The  sets  of  ao  values  were 
averaged  to  assign  a  single  ao  to  eaeh  loeal  piteh  angle 
bin,  and  time  interval  (/[,  ti).  Sehematieally,  these  steps  are 

(x.B)(t)  ^  {K,L}(t)  {ao.L}(t)  -  0) 

where  the  braees  indieate  a  set  eorresponding  to  the  18 
different  values  of  loeal  piteh  angle  (\. 

[20]  Flux  measurements  taken  during  eaeh  interval  (/i, 
^2)  were  identified,  and  log(/)  was  time  averaged  for  eaeh 
value  of  a  and  E  to  uniquely  speeify  j  as  a  funetion  of  ao, 
E,  and  7  =  (/|  +  Next,  sinee  the  data  were  far  too 
sparse  to  eover  all  values  of  ao,  the  flux  was  fit  to  a 
funetion  of  the  form  A  sin^'  ao  for  eaeh  MEA  energy 
ehannel  (or,  if  n  was  negative,  j  was  simply  averaged, 
equivalent  to  setting  n  =  0).  Interpolating  in  ao,  where  the 
data  were  suffieient,  yields  the  flux  values  shown  in  Figure  6. 
In  making  this  plot  of j\L,  /),  the  eonstant  value  J2  =  1 .78  x 
10  g(em/s)/?^  was  ehosen.  This  follows  Brautigam  and 

Albert  [2000],  who  were  performing  1-D  simulations  of 
y(L,  /)  at  eonstant  values  of  and  J2  and  determined  that 
this  value  of  J2  maximized  the  overlap  with  the  available 
data.  This  should  remain  roughly  true,  even  though  a 


different  magnetie  field  model  is  used  here,  and  allows 
for  at  least  rough  eorrespondenee  to  the  previous  work. 

[21]  The  data  eoverage  was  then  extended  by  extrap¬ 
olating  the  A  sin"  ao  fits  beyond  the  measured  values 
and  by  linear  interpolation  and  extrapolation  of  log  j  in 
log  E.  Finally,  for  arbitrary  times  /,  linear  interpolation 
of  log  j{t)  was  performed  at  fixed  (Aq,  L).  This 
allows  j  to  be  evaluated  at  any  grid  point  (a„  £),  L*)  at  any 
time.  Sehematieally, 

\{J}W  [{y}l(^i«^2)  [y(<>o)l(7)  — y(o/,£:,r) 

^j{ai.Ej,Lkj),  (8) 

where  the  braekets  refer  to  a  set  eorresponding  to  the  17 
different  energy  ehannels  of  MEA.  The  results  are  shown  in 
Figure  7,  for  the  same  fixed  value  of  J2.  (No  values  arc 
shown  below  L  =  4.2  for  0.42  MeV,  because  this  is  off  the 
computational  grid.)  These  fits,  interpolations,  and  extra¬ 
polations  of  the  available  data  are  regrettable  but  unavoid¬ 
able  if  values  are  to  be  determined  (or  assigned)  to  the  entire 
computational  grid. 

[22]  Onee  determined,  the  fluxes  were  converted  to  phase 
space  density  and  are  shown  in  Figure  8  at  fixed  values  of 
first  adiabatie  invariant  M  =  J\  (given  in  units  of  MeV/G). 
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Figure  6.  Electron  flux  as  measured  by  CRRES  MEA,  in  units  of  number  per  (cm~  s  sr  MeV),  as  well 
as  Dst  and  Kp. 


These  values  are  now  taken  to  represent  the  actual  data  and 
will  be  used  to  initialize  the  simulations,  to  drive  them  at  the 
boundaries,  and  for  comparison  to  the  results. 

4.  Simulations 

[23]  Data  from  CRRES  MEA  were  processed  as  just 
described,  for  the  interval  8-18  October  1990  (day  of  year 
281-291).  Dst  for  that  interval  is  shown  in  Figure  6  and 
indicates  a  moderate  geomagnetic  storm  beginning  during 
9  October  (day  282  of  1990).  As  discussed  in  detail  by 
Brautigam  and  Albert  [2000],  a  storm  sudden  commence¬ 
ment  at  1315  UT  (time  282.54)  was  accompanied  by 
strong  flux  decrease  at  L  =  5  for  £"  =  0.42-1.47  MeV. 
This  was  followed  by  an  injection  of  several  hundred  keV 
electrons  at  L  >  6  an  hour  into  the  recovery  phase  (time 
283.37).  Within  the  next  5  h  (by  time  283.58),  ~100  keV 
electron  flux  at  L  =  3-6  had  greatly  increased,  while 

1  MeV  flux  increased  much  more  gradually.  This  trend  of 
increase  near  L  =  5  was  interrupted  by  notable  dips  around 


t  =  285.0  and  t  =  288.5.  This  behavior  is  reproduced  in  the 
interpolated  and  extrapolated  data  of  Figure  7. 

[24]  The  corresponding  values  shown  in  Figure  8  were 
used  not  only  to  initialize  the  simulation  but  also  for  time- 
dependent  boundary  values  on  all  six  planar  faces  of  the 
3-D  simulation  domain.  Thus,  both  radial  diffusion  and 
local  heating  were  supplied  with  realistic,  dynamic  “seed 
populations”  from  which  to  generate  flux  at  relativistic 
energies.  Grid  resolution  was  43  values  of  Q\,  25  points 
in  Qi^  and  27  points  in  L,  covering  3.55  <  L  <  6.15  with 
AL  =  0.10,  Aao  ^  2°,  and  Ej+dEj  1.1  (although  only  in 
L  was  the  spacing  constant,  as  discussed).  For  simplicity,  a 
fully  explicit  finite  differencing  scheme  was  used,  which 
limited  the  time  step  by  the  Courant-Friedrichs-Lewy 
condition  for  linear  stability,  namely, 
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Figure  7,  CRRES  MEA  electron  flux,  interpolated  and  extrapolated  to  cover  the  computational  grid. 


7  of  21 


=  1000  MeV/G  M=500  MeV/G  M=200  MeV/G  M=100  MeV/G 


A09214 


ALBERT  ET  AL.;  THREE-DIMENSIONAL  DIFFUSION  SIMULATION 


A09214 


282  284  286  288  290 

t  (day  of  1990) 

Figure  8.  CRRES  MEA  electron  flux,  interpolated,  extrapolated,  and  converted  to  phase  space  density / 
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The  actual  time  step  was  taken  to  be  0.5 Atcri  ^nd  was 
reevaluated  whenever  Kp  changed.  This  resulted  in  Ar  1 6  s 
for  Kp  <  2,  At  ^  4.3  s  for  2  <  AT/?  <  4,  and  At  ^  0.77  s 
for  AT/?  >  4.  The  entire  9.5  day  simulation  ran  in  about  50  min 
on  a  standard  PC. 

4.1.  Overview 

[25]  As  mentioned,  the  value  J2  =  1  -  78  x  10  *^  g(cm/s)/?e 
was  used  for  the  comparisons.  First,  the  wave-induced 
pitch  angle  and  energy  diffusion  were  omitted,  leaving  just 
radial  diffusion.  Results  are  shown  in  Figure  9  and 
qualitatively  reproduce  the  results  of  Brautigam  and  Albert 
[2000].  In  particular,  fairly  good  agreement  with  the  meas¬ 
urements  (shown  in  Figure  8)  is  found  for  A/  =  100  and 
A/  =  200  MeV/G,  although  the  dropouts  around  t  =  283 
arc  too  weak,  while  the  increases  for  M  =  500  and  M  = 
1000  MeV/G  starting  around  t  =  286  are  far  too  weak 
and  transient.  As  noted  by  Brautigam  and  Albert  [2000], 
the  results  are  essentially  driven  by  transport  of  the  time- 
dependent  values  at  the  outer  radial  boundary,  L  =  6.15. 
As  found  previously,  this  is  sufficient  to  account  for  the 
observed  increases  at  lower  L  for  M  <  200  but  evidently 
not  for  A/  >  500. 

[26]  Next,  a  simulation  was  done  with  diffusion  in  (ao,  E) 
but  omitting  radial  diffusion.  As  shown  in  Figure  10,  this 
leads  to  large,  widespread,  sustained  increase  in  phase  space 
density  for  A/  >  200,  far  larger,  in  fact,  than  seen  in  the  data, 
especially  at  L  >  4.5.  Finally,  allowing  diffusion  in  ao,  E, 
and  L  to  operate  leads  to  intermediate  values,  as  seen  in 
Figure  11.  Varotsou  et  a  I  [2008]  also  found  this  ordering  of 
phase  space  density  values  when  including  chorus  and 
radial  diffusion  separately  and  together.  Figure  1 1  also 
shows  that  these  intermediate  values  also  match  the  data 
of  Figure  8  fairly  well. 

4.2.  Detailed  Evolution 

[27]  A  more  detailed  comparison  is  shown  in  Figure  12, 
which  shows  j[t)  from  the  three  simulations  at  A/  =  200 
and  M  -  1000  MeV/G,  at  several  fixed  values  of  L. 
Results  for  A/  =  100  are  similar  to  those  for  M  =  200, 
and  results  for  M  =  500  are  similar  to  those  for  M  =  1000. 
(The  corresponding  values  of  E  and  (Vq  can  be  determined 
from  the  top  row  of  Figure  5.) 

[28]  In  all  cases,  the  dropout  around  t  -  283  is  not  fully 
captured,  especially  at  low  L  (hence  high  £),  presumably 
because  of  wave-induced  precipitation  not  represented  in 
the  simulations,  such  as  by  hiss  (in  the  plasmasphere  and  in 
plumes),  electromagnetic  ion  cyclotron  waves,  or  possibly 
by  fast  magnetosonic  waves  [Li  et  ai,  2007;  Horne  et  ai, 
2007;  Albert,  2008].  A  check  verifies  that  /  at  A/  =  1000 
does  decrease  rapidly  near  the  dropout  for  larger  values  of 
J2  (smaller  values  of  rvo)  in  response  to  lower  values  at  the 
corresponding  grid  boundary,  but  evidently  the  pitch  angle 
diffusion  rates  are  too  low  for  the  values  shown  to  change 
much  before  the  boundary  conditions  recover.  Of  course, 
this  is  subject  to  limitations  in  deriving  the  boundary 
conditions  from  fits  to  the  limited  data. 

[29]  On  the  other  hand,  the  increases  are  captured  rather 
well  by  the  combination  of  radial  and  chorus  diffusion, 
which  works  better  than  either  mechanism  acting  alone.  The 
largest  discrepancies  arc  for  M  =  1000  at  A  =  4.55  and 


especially  2X  L  -  4.05,  where  the  small  values  (at  the 
dropout)  are  far  too  large,  and  the  large  values  (late  in  the 
simulation)  are  too  small  by  about  a  factor  of  2  or  3.  For 
L  >  4.55,  the  chorus  and  results  usually  lie  below  the 
chorus-only  values  and  above  the  D/ /.-only  values.  This  was 
evident  in  the  2-D  plots.  Thus,  chorus  seems  to  act  as  the 
source  of  phase  space  density,  while  radial  diffusion  acts 
mostly  to  transport  it  away.  However,  at  A  =  4.05  the  values 
from  combined  chorus  and  radial  diffusion  are  higher  than 
from  either  alone,  which  implies  net  radial  diffusion  into,  not 
away  from,  L  =  4.05. 

[30]  It  is  reasonable  to  question  the  development  of 

agreement  at  late  time  (/  290)  from  substantial  disagree¬ 

ment  at  earlier  times  (/  «  283).  Therefore,  the  simulations 
were  repeated  starting  at  r  =  283.4,  when  the  measured 
fluxes  arc  near  minimum.  The  results  are  shown  in  Figure  13 
and  are  generally  seen  to  revert  to  the  same  values  as  the 
previous  run  after  a  day  or  two.  This  indicates  that  fluxes  arc 
determined  by  transport  of  the  time-dependent  sources  at  the 
boundaries  more  than  by  existing  interior  values.  Varotsou 
et  al.  [2008]  also  found  that  large  differences  in  initial 
conditions  could  lead  to  relatively  similar  states  after  about 
a  day.  This  is  not  surprising,  since  the  time  scales 
associated  with  the  diffusion  coefficients  are  ~1  day,  as 
shown  in  Figure  3.  An  exception  occurs  for  M  =  1000  at 
L  -  4.05  for  the  chorus-only  run,  which  yields/ almost  as 
low  as  the  Du-only  run.  Here  it  is  very  evident  that 
chorus  and  Z)//.  do  not  compete,  but  instead  cooperate,  to 
produce  the  recovery  of  phase  space  density. 

4.3.  Sensitivity’  to  Diffusion  Coefficients 

[31]  As  a  sensitivity  test,  the  radial  and  chorus  diffusion 
coefficients  were  both  included  but  were  doubled  and 
halved,  separately  and  together.  Starting  the  simulation  at 
t  =  281.5,  the  effects  were  small  for  M  =  200  MeV/G 
(especially  at  large  L)  and  substantial  for  A/=  1000  MeV/G, 
as  seen  in  Figure  14.  Relative  to  the  “standard”  run,  shown 
again  as  the  black  curves,  increasing  (solid  blue  curve) 
led  to  a  more  realistic,  though  lagged,  dropout  around 
t  =  283.4  and  lower  recovered  values  of  /  at  late  times. 
Decreasing  Dchoms  (dashed  red  curve)  had  a  similar  effect, 
while  increasing  Dchorus  (solid  red  curve)  or  decreasing 
Di^i  (dashed  blue  curve)  tended  to  have  the  opposite 
effect,  leading  to  larger  /  both  at  the  dropout  and  later. 
Doubling  or  halving  the  strength  of  both  processes  at  the 
same  time  (solid  and  dashed  green  curves,  respectively) 
tended  to  produce  smaller  changes,  suggesting  that  chorus 
and  radial  diffusion  compete  in  determining  /  However, 
this  interpretation  is  not  consistent  with  the  runs  starting  at 
t  =  283.4,  shown  in  Figure  15,  especially  for  low  L  and 
large  A/.  Decreasing  the  strength  of  either  process,  or  both, 
led  to  considerably  lower  /  again  indicating  that  here 
chorus  and  radial  diffusion  reinforce  each  other  in  produc¬ 
ing  the  reeovery  of  f  Also  note  that  in  both  Figures  14 
and  15,  the  runs  with  doubled  D^homs  produced  excellent 
agreement  with  the  measured  values,  except  for  the  shal¬ 
lowness  of  the  dropout  obtained  in  Figure  14. 

4.4.  Radial  Profiles 

[32]  Figures  16  and  17  show  the  simulation  results  as 
radial  profiles,  f{L)  at  fixed  M  and  J2.  Figure  16  shows 
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Figure  9.  Phase  space  density,  simulated  with  radial  diffusion  only. 
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Figure  10.  Phase  space  density,  simulated  with  chorus  diffusion  only. 
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Figure  11.  Phase  space  density,  simulated  with  both  radial  diffusion  and  chorus  diffusion. 
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Figure  12.  Phase  space  density,  in  units  of  s*^  km  as  determined  from  CRRES  data  (black  diamonds), 
and  simulations  starting  at  t  -  28 1.5,  with  only  (blue  curves),  with  chorus  only  (red  curves),  and  with 
both  (black  curves). 


13  of  21 


=4.05  L=4.55  L=5.05  L=5.55 


A09214 


ALBERT  ET  AL.:  THREE-DIMENSIONAL  DIFFUSION  SIMULATION 


A09214 


10 


-8 


10 


-9 


10 


-10 


10 


-10 


282 


284  286  288  290 

t  (doy  of  1 990) 


Figure  13.  Phase  space  density,  in  units  of  km“^,  as  determined  from  CRRES  data  (black  diamonds), 
and  simulations  starting  at  /  =  283.4,  with  only  (blue  curves),  with  chorus  only  (red  curves),  and  with 
both  (black  curves). 
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Figure  14.  Phase  space  density  as  determined  from  CRRES  data  (black  diamonds)  and  simulations 
starting  at  /  =  281.5,  with  both  and  chorus  (black  curves),  with  Du  doubled  (solid  blue  curves)  and 
halved  (dash-dotted  blue  curves),  /^chorus  doubled  (solid  red  curves)  and  halved  (dash-dotted  red  curves), 
and  with  both  doubled  (solid  green  curves)  and  halved  (dash-dotted  green  curves). 
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Figure  15.  Phase  space  density  as  determined  from  CRRES  data  (black  diamonds)  and  simulations  as 
in  Figure  14,  starting  at  t  =  283.4. 
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Figure  16.  Snapshots  of  phase  space  density  versus  L  at  various  values  of  t  for  M  =  1000  MeV/G, 
simulated  with  Dll  only,  Dghorus  only,  and  both.  CRRES  data  at  the  end  of  the  simulation  is  shown  as  red 
diamonds  in  all  three  plots. 


snapshots  at  several  values  of  t  for  M  =  1000  MeV/G, 
from  the  start  of  the  simulations  at  t  =  281.5  to  the  end 
near  /  =  291,  as  well  as  the  CRRES  data  at  the  ending 
time.  The  initial  profile  inereases  essentially  monotonieally 
with  inereasing  L.  The  simulation  with  only  radial  diffu¬ 
sion  develops  internal  peaks,  eaused  by  the  varying 
boundary  eonditions  at  the  maximum  and  minimum  L, 
but  ends  up  again  monotonie  and  maximum  at  the  outer¬ 
most  L,  in  qualitative  as  well  as  quantitative  disagreement 
with  the  data.  The  chorus-only  simulation  produces  pro¬ 
files  that  remain  monotonie  or  nearly  so,  as  well  as 
beeoming  too  large.  This  is  contrary  to  what  might  be 
expected  from  a  localized  internal  souree  [Green  and 
Kivelson,  2004].  The  eombination  of  ehorus  and  radial 
diffusion  produecs  robust  internal  peaks  around  L  -  4.5 
that  resemble  the  data,  though  they  are  a  bit  too  low  as 
noted  in  seetion  4.2. 

[33]  Figure  17  shows  analogous  snapshots  of  /(L)  for 
several  values  of  A/.  The  eurves  in  Figure  17a,  based  on  the 
fit  CRRES  data,  start  out  monotonieally  increasing  with 
increasing  L  but  quickly  develop  internal  peaks  at  all  M 
which,  for  M  =  500  and  M  -  1000  MeV/G,  persist 
throughout  the  time  interval.  The  simulations  with  Dll^ 
with  or  without  chorus,  reproduce  this  behavior  for  A/  =  100 
and  M  =  200  MeV/G,  but  without  chorus  the  simulated 
peaks  at  A/  =  500  and  M  -  1 000  are  too  short  lived.  The 
chorus-only  run  never  displays  internal  peaks.  Only  the  run 
combining  chorus  with  Dll  shows  qualitative  agreement 


with  the  data  throughout  the  run,  at  both  low  and  high 
values  of  A/. 

4.5.  Effect  of  Cross  Terms 
[34]  Finally,  we  briefly  consider  the  effect  of  omitting  the 
cross-diffusion  terms,  involving  D^^^p.  The  placement  of 
grid  points  in  (ao,  £)  was  unchanged,  but  the  coefficient 
Da^p  was  artificially  set  to  0,  and  the  numerical  procedures 
(such  as  tracing  constant-02  curves,  which  become  con- 
stant-E  curves)  were  carried  out  as  before.  Figure  1 8  shows 
the  ratio  ofy(ao,  E)  with  and  without  cross  diffusion  at  /  = 

285.5  and  L  ~  4.05  for  the  chorus-only  run  starting  t  = 
283.4.  As  expected  from  previous  2-D  studies  with  similar 
chorus  models  [Albert  and  Young,  2005;  Tao  et  ai,  2008, 
2009],  the  effect  is  concentrated  above  E  -  1  MeV  and  low 
values  of  ao,  with  a  peak  ratio  of  about  50,  and  is  modest 
elsewhere.  Figure  19  shows  how  this  ratio,  evaluated  at 
E  =  1.8  MeV  and  Oq  =  25®,  behaves  with  time  (dash-dotted 
red  curve).  It  is  initially  1  (since  the  initial  conditions  are 
identical  regardless  of  cross  terms),  but  quickly  grows  to  a 
persistent  level  of  about  10.  This  case,  which  simulates  the 
production  of  electrons  of  energy  ~  1  MeV  from  a  depleted 
state,  is  representative  of  how  2-D  simulations  of  storm  time 
chorus  heating  are  usually  performed.  However,  for  the 
chorus-only  run  starting  at  /  =  281.5  (red  solid  curve),  the 
ratio  only  grows  to  about  3.  More  significantly,  the  runs  with 
both  chorus  and  radial  diffusion  also  show  a  ratio  of  about  3, 
starting  from  either  t  -  283.4  (dashed  black  curve)  or  t  = 

281.5  (solid  black  curve).  Thus,  at  least  for  the  diffusion 
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Figure  17.  Snapshots  of  phase  space  density  versus  L  at  various  values  of  t,  determined  from 
(a)  CRRES  data,  and  (b)  simulated  with  Du.  only,  (c)  Ochoms  only,  and  (d)  both. 


models  used  here,  the  substantial  effect  of  the  cross  terms  in 
2-D  simulations  is  considerably  reduced  in  the  presence  of 
radial  diffusion. 

5.  Summary  and  Discussion 

[35]  In  summary,  we  have  performed  a  3-D  simulation  of 
the  9  October  1990  magnetic  storm  accounting  for  radial 


diffusion  and  chorus  wave-induced  diffusion  in  pitch  angle 
and  energy,  including  cross  terms.  Grid  points  were  aligned 
in  a  natural  way  for  L  diffusion  (at  constant  M  and  JiX  while 
at  each  L  a  grid  in  (Qi,  Q2)  was  constructed  for  the  chorus 
diffusion.  In  the  spirit  of  Brautigam  and  Albert  [2000], 
CRRES  MEA  particle  data  have  been  used  to  obtain 
realistic,  time-dependent  boundary  conditions,  which  are 
transported  throughout  the  computational  domain  using 
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Figure  18.  Ratio  of  simulated  /without  cross  ditTiision  to /with  cross  diffusion,  at  Z.  =  4.05,  /  =  285.5, 
for  the  chorus-only  run  starting  at  t  =  283.4. 

activity-dependent  diffusion  coefficients.  The  resulting 
agreement  with  data  is  reasonably  good,  especially  when 
chorus  diffusion  is  increased  by  a  factor  of  2,  except  for  the 
depth  of  the  main  phase  dropout.  However,  this  effect  of 
this  deficiency  is  short  lived,  since  runs  initialized  before 
and  immediately  after  the  dropout  quickly  converge  to  each 


other.  (The  beneficial  factor  of  2  for  should  not  be 

taken  too  literally,  given  the  uncertainties  in  the  wave 
models.) 

[36]  Experimentation  shows  that  the  3-D  combination  of 
Dll  snd  chorus  is  essential  and  that  either  process  alone 
does  not  give  even  rough  quantitative  agreement  with  the 


Figure  19.  Ratio  of  simulated / without  cross  diffusion  to / with  cross  diftusion,  at  L  =  4.05,  with  chorus 
only  (red)  curves  and  with  both  Du  and  chorus  (black  curves),  starting  at  /=  28 1 .5  (solid  curves)  and 
at  r  =  283.4  (dash-dotted  curves). 
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persistent  internal  peaks  seen  in  the  data.  Such  phase  space 
density  peaks  are  a  common  feature  seen  during  magnetic 
storms  [e.g.,  Green  and  Kivelson,  2004;  lies  et  aL,  2006; 
Chen  et  ai,  2007],  Furthermore,  the  combination  of  the  two 
processes  is  complex,  since  chorus  can  cause  either  increase 
phase  space  density  through  energy  diffusion  or  decrease  it 
by  pitch  angle  diffusion,  and  radial  diffusion  can  act  to 
either  increase  or  decrease  /depending  on  gradients.  Thus, 
simple  interpretations  based  on  “competition,”  or  indeed 
“cooperation,”  can  be  misleading.  The  three-dimensional 
simulations  presented  here  support  the  paradigm  of  inward 
radial  diffusion  of  lower-energy  “seed”  electrons  which  are 
energized  by  chorus  waves,  and  then  radially  diffused  both 
inward  and  outward,  resulting  in  the  observed  internal  peaks 
[Horne,  2007]. 

[37]  As  indicated,  the  reasonable  success  in  reproducing 
the  CRRES  data  for  this  storm  depended  on  having  bound¬ 
ary  values  on  all  six  of  the  grid  boundaries,  which  raises  the 
question  of  practicality  for  space  weather  forecasting.  While 
it  is  not  unreasonable  to  presume  the  availability  of  data  at 
all  needed  values  of  (ao,  E)  at  L^nm  and  one  cannot 
count  on  having  a  time  series  measurements  at,  say,  all 
(ao,  L)  at  fixed  Emin-  However,  these  may  be  supplied  by  a 
ring  current  code,  which  almost  by  definition  aims  to 
simulate  particles  up  to  lower  radiation  belt  energies. 

[38]  These  results  support  the  effectiveness  of  simulating 
chorus-electron  interactions  by  quasi-1  inear  diffusion,  de¬ 
spite  the  increasingly  appreciated  nonlinear  nature  of  chorus 
waves.  More  work  is  therefore  needed  not  only  to  develop 
quasi-1  inear  modeling  but  to  understand  why  it  seems  to 
work  as  well  as  it  does. 
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